Self-organized criticality in stick-slip models with periodic boundaries 
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A spring-block model governed by threshold dynamics and driven by temporally increasing spring 
constants is investigated. Due to its novel multiplicative driving, criticality occurs even with periodic 
boundary conditions via a mechanism distinct from that of previous models. This mechanism is 
dictated by a coarsening process. The results show a high degree of universality. The observed 
behavior should be relevant to a class of systems approaching equilibrium via a punctuated threshold 
dynamics. 

PACS numbers: 64.60.Lx, 64.60.-i, 05.70.Jk 



Out-of-equilibrium driven systems with threshold dy- 
namics exhibit a rich phenomenology, from synchronized 
behavior to self-organized criticality (SOC) |^-^. 
SOC refers to the spontaneous organization towards a 
kind of dynamical critical steady state. Threshold out-of- 
equilibrium dynamics encompasses many systems, such 
as neural networks, solid friction, rupture with healing, 
earthquakes and avalanches. It is now understood qual- 
itatively that there is a class of models exhibiting SOC 
as a result of their tendency to synchronize . This ten- 
dency is however frustrated by constraints such as open 
boundary conditions and quenched disorder Q which 
leads to a dynamical regime at the edge of synchroniza- 
tion, the SOC state. Another class, so-called extremal 
models, are understood to exhibit SOC due to the com- 
petition between local strengthening and weakening due 
to interactions ||]. In a third class of models, SOC re- 
sults from the tuning of the order parameter of a system 
exhibiting a genuine critical point to a vanishingly small, 
but positive value, thus ensuring that the corresponding 
control parameter lies exactly at its critical value for the 
underlying depinning transition |^ . The issue is further- 
more complicated by the fact that a notable fraction of 
numerical and experimental works |^^|-|l2| claiming the 
observation of SOC from the measurements of powerlaw 
distributions rely on the slow sweeping of a control pa- 
rameter towards a critical point 
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The purpose of this letter is to present a variation of 
spring-block models using a novel form of driving by mul- 
tiplicative loading. The surprising result is that, when the 
dynamics is conservative and for periodic boundary con- 
ditions (PBC), the system self-organizes into a critical 
state with long-range order and power-law distributions. 
This is in contrast to all previous stick-slip models P,p|p^ 
that need either open boundary conditions and/or dissi- 
pation to exhibit SOC (in order for synchronization not 



to occur). None of the four above mechanisms seem at 
work here. We conjecture that this new class appears as 
a result of the multiplicative nature of the driving, known 
in simpler models to lead to stationary dynamical states 
with power law distributions u& . Our model exhibits the 
dual property of approaching (without ever reaching) an 
equilibrium state in terms of block positions while at the 
same time being characterized by a steady state in the 
force variables. To our knowledge our model is the first 
of its kind to show SOC in approach to an equilibrium 
state, and as such suggests a new class of experimental 
systems which could exhibit SOC states. 

Model: We consider a two-dimensional spring-block 
model consisting of an array of blocks interconnected 
among nearest neighbors by coil springs. The springs 
have the same spring constant K and relaxed spring 
length I. Initially, the array is stretched to a lattice spac- 
ing a > I and placed on a frictional substrate which is 
characterized by a static threshold Fs for slipping. Dis- 
order is introduced in the form of random displacements 
(x, y) of the blocks about the coordinates {ia, ja) on a 
square lattice, where —A < x,y < A and i,j — 1, .... L. 
The force between two neighboring blocks at r and r' 
is given by the Hooke's law (|r — r'\ — l)K. Since we 
are interested in dynamics primarily governed by tensile 
stresses, this nonlinear dependence on the coordinates 
leads to unnecessary complications in the algorithm. To 
simplify and compare with similar models, we expand 
the expression to first order in {x,y) to obtain the force 
components on a block in the bulk at {i,j) 

^ {xi+i.j + Xi^ij - 2xi^j)K 

+ {X'i,j+i + Xi^j^i - 2xi^j)sK 
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where s = l — l/a>Ois the initial strain. It is impor- 
tant to stress that the terms proportional to s leads to 
anisotropic couplings to nearest neighbors in the trans- 
verse direction. The coupling and the SOC state disap- 
pears for s = or in ID-chains. 

Since the forces are linear in displacements, it is possi- 
ble to invert (|l|) and formulate the model solely in terms 
of force variables, as in Q]. Starting with a stable con- 
figuration with net force F = \/ F^"^ + < Fg for all 
the blocks, wc drive the system by gradually increasing 
K 18 1 until one of the blocks becomes unstable, i.e., K 
is increased to KFs/Fniax during loading, where Fmax 
denotes the spatial maximum of F in the stable con- 
figuration. As in Q, the block is assumed to slip to its 
equilibrium position defined by -F = 0, ignoring overshot: 
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where a = 1/(2 + 2s). This locally conserves the force 
components, thus the conservative nature of the model. 
The resulting modification of the stress environment may 
trigger further slips in neighboring blocks, and hence an 
avalanche, until F < Fg is restored for all blocks. Then 
K is increased again and the slip process continues. 

If Fg was zero, the only stable (minimal energy) con- 
figuration would be that the blocks were exactly at the 
nodes of a perfect square lattice of mesh size a. The 
non-vanishing friction thus creates a large ensemble of 
coexisting metastable states which is responsible for the 
nontrivial dynamics. For Fg and s ^ 0, the toppling 
rules (||) do not put the blocks in their minimum energy 
configuration due to the couplings to their four neighbors. 
This ensures that a block will go on becoming unstable 
ad infinitum as long as K is increased indefinitely. 

The multiplicative loading is motivated by the stiff- 
ening of an overlayer caused by desiccation ]T9[|, origi- 
nally used to study cracks j2^. It differs from the usual 
additive loading in sandpiles |^-^ and stick-slip mod- 
els fl^Jl^ where physically the driving force arises from 
the dropping of grains onto a pile or from the steady 
relative motion of two frictional surfaces. Those sys- 
tems are known to exhibit self-organized criticality in 
non-equilibrium steady states. In contrast, our system 
approaches an equilibrium instead of a genuine steady 
state. 

Without loss of generality, we hereafter set a = 1 = . 
Of the remaining dimensionless parameters {s,A,K}, s 
determines the equilibrium length scale and the dynamics 
through (|^), A characterizes the initial disorders but is 
irrelevant for the equilibrium state, and K defines the 
"time" t = K. 

Results: We have simulated the model for < s < 1, 
A = 0.02,0.18,0.4, and system sizes 20 < L < 300. We 
are mainly interested in the possibility of SOC with PBC 
as this most differentiates us from previous works. Unless 



stated otherwise, PBC is assumed hereafter. To investi- 
gate the effects of spatial inhomogeneities, we also use 
free boundary conditions (FBC) with no block beyond 
the edges and corners, and cylindrical boundary condi- 
tions (CBC) with one pair of parallel edges periodic and 
the other pair free. All these boundary conditions respect 
the conservation law of the force component. 

The evolution of the system is monitored by the vari- 
ation of the stress field a{r,t), approximated by the av- 
eraged tension of the four springs attached to a block. 
Using its Fourier transform tr, we compute the struc- 
ture factor 5(fc,t) = {\a{k,t)\^) / - L'^Sj: ^{d{k,t)f 

and the circular average Sci-[(k,t) — X]fce{|fc|=fc} '^(^' 
where the overline and the angular brackets mean a 
spatial and an ensemble average, respectively. Then 

■R(^) = Efe k'^Sc„{k, t)l Y,k '5cir(fc, t)] gives a simple 
measure of the characteristic length for the stress field. 

The system evolves under a conservative dynamics 
from large disorders to small disorders when the blocks 
converge onto a perfect lattice. This is analogous to 
spinodal decomposition ||2l| and suggests a coarsening 
in the stress field. We indeed find a power-law growth 
R(t) ~ t"^, where = 0.33 ± 0.01 (cf. Q) for aU s and 
A, as illustrated in Fig. 1. This universal behavior and 
the value of (/> agree with spinodal decomposition. 

Now we present the evidence showing the system in the 
long-time limit is stationary and critical in the variables 
relevant to the dynamics. Fig. 1 shows the first evidence 
of stationarity in i?(t), where it reaches a plateau after 
transcient. While the saturation rate depends on s and 
A, the value of the plateau i?o depends only on s due 
to sL being the length the system has to contract to 
reach a stress-free state. But this cannot be achieved 
due to PBC, so that the blocks wind up in a frustrated 
state correlated over this distance. For fixed i, we ver- 
ify that Ro{s, L) increases linearly in s, except very close 
to s = or 1. The plateau extends up to ten decades 
in t, until reaching the limit of numerical accuracy (i.e., 
10^^^ in F, using double precision). Further evidence 
comes from F(t)/Fs which measures the effective "dis- 
tance" from the instability limit for the system. The 
stationary fluctuations of this ratio about a finite con- 
stant, as shown in Fig. 2, is an important characteristic 
of a dynamical steady state in SOC models 

Next, we show that the system is critical. Firstly, for 
fixed s, we find Rq{s,L) cx L, implying long-range cor- 
relations in the stress field. More importantly, R satis- 
fies finite-size scaling R{t,L) = R{t / L'^ / (see Fig. 3), 
with the asymptotic behavior R{x 0) = const and 
R{x — > oo) (X a;^*^. This is a clear signature of criticality. 

Secondly, the avalanches are characterized by power- 
law distributions, another hallmark of criticality. From 
(|l|), the force drop in one block-slip is given by = 
Ku/a PJT^], where u denotes the slip distance. Since 
^shp > p^ y^ diminishes as 1 /K when the blocks gradu- 
ally converge to a regular lattice. An avalanche consists 
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of a correlated sequence of S block-slips. Analogous to 
seismology the "seismic moment" M and the "radi- 
ated seismic energy" E can be computed: 
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These can be measured experimentally, as is precisely 
done for earthquakes. In our context, E equals the en- 
ergy dissipated by friction. We find that the distribu- 
tion of 5* approaches a power law P{S) ~ with 
B = 0.15 ± 0.05 (see Fig. 4), when R{t,L) Ro{L) 
for t > L^/"^. Consequently, both distributions of M 
and the scaled energy EK are stationary and follow the 
same power law. From these results, we conclude that 
our system approaches in a punctuated manner without 
ever reaching its final equilibrium, and is driven into a 
marginally stable state which is critical in physical vari- 
ables. 

The approach to thermodynamic equilibrium via spin- 
odal decomposition |2l] differs from our case in that 
its dynamics is governed by thermal fluctuations, not a 
threshold. On the other hand, it is well known that all 
previous SOC models constructed with threshold 

dynamics do not exhibit SOC with PBC, because the 
periodic boundaries forbid an outflux of the conserved 
quantity (the total number of grains or the total force 
on the array) to balance its external influx. The key 
difference with our model is that the variables F^ and 
Fy which are redistributed conservatively can be of ei- 
ther sign, with the total force remaining at zero during 
both block slips and loading, so that F which determines 
topplings (i.e., slips) is neither monotonic nor conserved. 
This separation of the conservation law from the toppling 
condition allows for criticality despite PBC. Another way 
to rationalize our results is to notice that, as mentioned 
above, the mechanisms at work to produce powerlaws 
in the presence of multiplicative noise (amplification by 
multiplication followed by reinjection) Pq | might also be 
relevant here. 

We also test the robustness of the critical properties 
against boundary effects. We find the same exponent B 
within numerical uncertainty for PBC, FBC, and CBC. 
This is remarkable in view of their different equilibrium 
states and paths of approach, and the usual sensitive de- 
pendence of SOC systems on boundary conditions. 

Our model depends significantly on only one factor: 
the conservation. To show this, we introduce non- 
conservation by adding a term —{x,y)KK to {Fx,Fy) 
in (|l|), which may represent harmonic couplings of each 
block with another surface, as in earthquake models 
The parameter k, then quantifies the level of non- 
conservation, with a in (|) replaced by l/{2 + 2s + n) 
For not too smaU k(> 0.1),F(i)/F^ does not 
show any stationary regime, and P{S) is exponential for 



all choices of sampling intervals in t. Consequently, the 
SOC state for k = is lost except at most for a small 
interval near k, — 0. 

Conclusion: It has been suggested that SOC in sandpile 
and spring-block models arises from a de-sychronization 
mechanism that is initiated by inhomogeneities from a 
free boundary |^ . We have studied a system which shows 
SOC without such inhomogeneities. This reveals a differ- 
ent mechanism whereby correlations (or "self-organized" 
regions) gradually build up in the bulk via a coarsening 
process. The completion time r « L^/*^ (cf. Fig. 3) is 
characterized by the coarsening exponent (j), which should 
be compared with the invasion time in previous models. 
The associated power-law exponent is extremely robust 
(universal) with respect to the initial disorder, the initial 
strain, the size of the system and the type of boundary 
conditions. While exhibiting a well-defined steady state 
in the characteristic length and the average force, the 
model has also a transient nature when viewed in the slip 
distance toward the equilibrium and in the sweeping of 
the spring coefficient K. This teaches us that experimen- 
tal systems that appear transiently driven might in fact 
be stationary in the variable relevant to the dynamics, 
especially when converging to a fundamental equilibrium 
state. Search schemes and optimization techniques using 
the sweeping of a control parameter such as in simulated 
annealing to get access to the fundamental state or to the 
optimal solution might exhibit this kind of phenomenon 
in which the relaxation is characterized by a wide distri- 
bution of jumps. 
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FIG. 1. Characteristic length R{t) vs time t = K for 
different sets of s and A, showing the universal 1/3 power law 
and plateau. L = 40 with PBC. 
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FIG. 2. Spatial averaged force on a block normalized by 
the threshold vs time, showing approaches to stationarity. 
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FIG. 3. Finite-size scaling plot of the characteristic length 
in the stress field, R{t, L). s = 0.5, A = 0.18 with PBC. 
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FIG. 4. P{S) obtained from successive intervals in t to 
show the approach to power law as R{t) — » Ro. 
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